function error = posteriori_biharmonic(dofs,V,T,posT,d,c,qw,qx,qy,phi4)
% error = posteriori_biharmonic(dofs,V,T,posT,d,c,qw,qx,qy,phi4)
% now the estimator is only fit for biharmonic problems
error = zeros(posT,1);
for k = 1:posT
    %comput the local bending matrix
    V1=V(T(k,1),:);V2=V(T(k,2),:);V3=V(T(k,3),:); degree = d(k);
    J = (V2(2)-V3(2))*(V1(1)-V3(1)) - (V1(2)-V3(2))*(V2(1)-V3(1));
    a = (V2(2)-V3(2))/J; b = (V3(2)-V1(2))/J;
    p = (V3(1)-V2(1))/J; q = (V1(1)-V3(1))/J;
    LocK = phi4{degree,1}*(a^4 + p^4 + 2*a^2*p^2) + ...
           phi4{degree,2}*(4*a^3*b + 4*p^3*q + 4*a*b*p*p + 4*a*a*p*q) +...
           phi4{degree,3}*(6*a*a*b*b + 6*p*p*q*q + 8*a*b*p*q + 2*b*b*p*p + 2*a*a*q*q) + ...
           phi4{degree,4}*(4*a*b^3 + 4*p*q^3 + 4*a*b*q*q + 4*b*b*p*q) + ...
           phi4{degree,5}*(b^4 + q^4 + 2*b^2*q^2);
    loc_dof = get_tri_dof(dofs,k);
    c_loc = c(loc_dof);
    % find longest side's length
    sides = [(V(T(k,2),1)-V(T(k,3),1))^2 + (V(T(k,2),2)-V(T(k,3),2))^2; ...
        (V(T(k,1),1)-V(T(k,3),1))^2 + (V(T(k,1),2)-V(T(k,3),2))^2; ...
        (V(T(k,1),1)-V(T(k,2),1))^2 + (V(T(k,1),2)-V(T(k,2),2))^2];
    longest =  max(sides);
    error(k) = J*qw'*(LocK*c_loc).^2*longest^2;
end